{"cells":[{"cell_type":"markdown","metadata":{"collapsed":false,"id":"6C09B47761DE44A58B01A2763B92E0EC","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"# Lecture 6: Denominator and Posterior\n\n## Instructor： 胡传鹏（博士）[Dr. Hu Chuan-Peng]\n\n### 南京师范大学心理学院[School of Psychology, Nanjing Normal University]"},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"EC6D8489DD764ABDABA8E3C20F4E21BB","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"## Part 2: The Posterior – The Goal of Bayesian Inference"},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"F22731CFF6AC43268CB888EEF4046926","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"\n![Image Name](https://cdn.kesci.com/upload/image/rjj2r4t0rq.png?imageView2/0/w/640/h/640)\n\n\n后验概率分布 $p(\\theta|data)$ 可以用来估计我们感兴趣的参数，比如疾病发生率，地球海水覆盖率等。\n\n在本节课程中，我们将探究后验概率分布是如何形成的，并讨论后验分布在实践中如何帮助我们解决问题。"},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"48FA13521A8C4E5185552CF3CDA3F046","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"### 后验如何更新通过数据更新我们的先验\n\n在贝叶斯统计中，后验分布将我们**已有的信念**与**观察到的数据信息**结合起来，产生一个**更新后的信念**，代表我们对知识更新的状态。\n\n为了计算后验分布，我们必须选择一个**似然函数**(第四节课)。似然函数的关键在于建立参数与数据的联系。\n\n在更新后验之前，参数需要首先被赋以**先验**，表示我们在得到数据前的信念(第三节课)。\n\n然后将先验和似然通过贝叶斯公式结合——得到后验分布:\n$p(\\theta|data) \\sim p(data|\\theta) * p(\\theta)$ \n\n注意，因为分母 $p(data)$是一个定值，在上一节课中提到它的作用在于使得后验分布成为一个有效(valid)概率分布。\n这里为了方便讨论先验和似然对于后验的影响，我们先忽略分母部分。"},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"DBC3EC62B4F84D6781F2C1D77C362A26","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"**以探究地球海水覆盖率为例。**\n\n- 假设我们在获得数据以前的**先验**认为：地球海水覆盖率为0-1的可能是相等的。\n\n\n"},{"cell_type":"code","execution_count":1,"metadata":{"collapsed":false,"id":"1029A8E8C3E24B67BBA6530F53EC7661","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[],"source":"import numpy as np               # numpy 是专门用于数组运算的包\nimport scipy.stats as st         # 从scipy.stats里载入分布函数\nimport matplotlib.pyplot as plt  # matplotlib.pyplot 是专门用于画图的包\nimport seaborn as sns            # seaborn是专门用于绘制统计数据的包\nimport arviz as az               # arviz 是专门用于探索贝叶斯模型的包"},{"cell_type":"code","execution_count":2,"metadata":{"collapsed":false,"id":"CE23D1B7C31745E78D9C198F7B94F695","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"output_type":"execute_result","data":{"text/plain":"Text(0, 0.5, 'probability density')"},"metadata":{},"transient":{}},{"output_type":"display_data","data":{"text/plain":"<Figure size 432x288 with 1 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/CE23D1B7C31745E78D9C198F7B94F695/rjj8jf3ibs.png\">"},"metadata":{"needs_background":"light"},"transient":{}}],"source":"x = np.linspace(0, 1, 100)\nprior = st.uniform.pdf(x, 0, 1)/np.sum(st.uniform.pdf(x, 0, 1))\n\nplt.plot(x, prior , color = 'grey',ls = '--')\nplt.xlabel(\"Sea surface coverage\")\nplt.ylabel(\"probability density\")"},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"AEE11F899B3C4C388132A96F19163EC1","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"- 之后，假设我们获得数据：抛地球仪10次，其中7次海面朝上。\n- 在获得数据的同时，我们以二项分布函数作为**似然函数**。"},{"cell_type":"code","execution_count":3,"metadata":{"collapsed":false,"id":"2D02930B13F5412297F10C5EAE5B48B8","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"output_type":"execute_result","data":{"text/plain":"Text(0, 0.5, 'likelihood')"},"metadata":{},"transient":{}},{"output_type":"display_data","data":{"text/plain":"<Figure size 432x288 with 1 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/2D02930B13F5412297F10C5EAE5B48B8/rjj8jkabts.png\">"},"metadata":{"needs_background":"light"},"transient":{}}],"source":"likelihood = st.binom.pmf(7, 10, x)\n\nplt.plot(x, likelihood , color = 'grey')\nplt.xlabel(\"Sea surface coverage\")\nplt.ylabel(\"likelihood\")"},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"5B35465A553D47C18154A63D86A95B46","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"- 通过贝叶斯公式 $p(\\theta|data) \\propto p(data|\\theta) * p(\\theta)$ ，得到后验分布。"},{"cell_type":"code","execution_count":4,"metadata":{"collapsed":false,"id":"27DEF2653AFA410390D8D42019B41B89","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"output_type":"execute_result","data":{"text/plain":"[<matplotlib.lines.Line2D at 0x7f92ac97d2d0>]"},"metadata":{},"transient":{}},{"output_type":"display_data","data":{"text/plain":"<Figure size 432x288 with 1 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/27DEF2653AFA410390D8D42019B41B89/rjj8jnp1ku.png\">"},"metadata":{"needs_background":"light"},"transient":{}}],"source":"unstd_posterior = prior*likelihood\n\nplt.plot(x, unstd_posterior , color = 'black')"},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"E2FA632B507042208E87304AB9F12685","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"在同一幅图中对比先验，通过可视化方式来比较似然与后验分布。"},{"cell_type":"code","execution_count":5,"metadata":{"collapsed":false,"id":"FBFF83C2C16C40E1A533548C99A1FEDA","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"output_type":"execute_result","data":{"text/plain":"<matplotlib.legend.Legend at 0x7f92ac896e10>"},"metadata":{},"transient":{}},{"output_type":"display_data","data":{"text/plain":"<Figure size 432x288 with 1 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/FBFF83C2C16C40E1A533548C99A1FEDA/rjj8jy1k90.png\">"},"metadata":{"needs_background":"light"},"transient":{}}],"source":"unstd_posterior = prior*likelihood\nposterior = unstd_posterior/np.sum(unstd_posterior) # 归一化\n\nplt.plot(x, prior, color = 'grey',ls = '--',label=\"prior\")\nplt.plot(x, likelihood, color = 'blue',label=\"likelihood\")\nplt.plot(x, posterior, color = 'red',label=\"posterior\")\nplt.legend()"},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"A61E7DAB7586478A9D0E53A1DEC09EF7","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"可以看到，由于先验设置得非常平滑（无先验信息），因此我们得到的后验分布也相对平滑。\n\n- 假设我们再收集10次抛地球仪的数据，得到6次海水向上。\n- 此时，我们可以将之前的后验作为先验再次通过数据更新我们的信念。"},{"cell_type":"code","execution_count":6,"metadata":{"collapsed":false,"id":"AFA211F27D42476F82AEECEB0DD8A5ED","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"output_type":"execute_result","data":{"text/plain":"<matplotlib.legend.Legend at 0x7f92ac85d050>"},"metadata":{},"transient":{}},{"output_type":"display_data","data":{"text/plain":"<Figure size 432x288 with 1 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/AFA211F27D42476F82AEECEB0DD8A5ED/rjj8k4eoce.png\">"},"metadata":{"needs_background":"light"},"transient":{}}],"source":"likelihood = st.binom.pmf(6,10, x)                   # 新数据的likelihood\nunstd_posterior = posterior*likelihood               # 将上次的后验当作这次数据的先验\nposterior2 = unstd_posterior/np.sum(unstd_posterior) # 进行归一化\n\nplt.plot(x, posterior , color = 'grey',ls = '--',label=\"prior\")\nplt.plot(x, likelihood, color = 'blue',label=\"likelihood\")\nplt.plot(x, posterior2, color = 'red',label=\"posterior\")\nplt.legend()"},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"E2E0476C523047988FFE9918E02B594A","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"可以看到，后验分布 (红色) 在先验 (灰色) 的基础上向似然更加接近了。\n\n这里体现了贝叶斯方法的一个优势：数据顺序不变性 (data order invariance)。\n- 即相对于频率学派通过一次收集数据进行统计推断。\n- 贝叶斯方法可以不断的获取数据来更新自己的推断。\n- 并且这种推断的过程不受数据收集先后顺序的影响。比如，第一数据为7次海面朝上，第二次数据为6次海面朝上."},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"C43A71F02FCA4224889C240C90252D20","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"甚至我们可以在每次获得数据后，对后验分布进行更新：\n\n\n![Image Name](https://cdn.kesci.com/upload/image/rhqd6akbc6.gif?imageView2/0/w/500/h/500)"},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"E6B81BE3FA9D47D089CDFE6F4BB3A8A7","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"我们可以把刚才多次对于后验的更新合并为一次："},{"cell_type":"code","execution_count":7,"metadata":{"collapsed":false,"id":"9054F9C77A814EF587BB8294A1F63E07","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"output_type":"execute_result","data":{"text/plain":"<matplotlib.legend.Legend at 0x7f92ac8222d0>"},"metadata":{},"transient":{}},{"output_type":"display_data","data":{"text/plain":"<Figure size 432x288 with 1 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/9054F9C77A814EF587BB8294A1F63E07/rjj8kbt9ei.png\">"},"metadata":{"needs_background":"light"},"transient":{}}],"source":"# 两次数据合并后，20次数据中有13次海水向上\nlikelihood = st.binom.pmf(13, 20, x)\nunstd_posterior = prior*likelihood\nposterior3 = unstd_posterior/np.sum(unstd_posterior)\n\nplt.plot(x, prior , color = 'grey',ls = '--',label=\"prior\")\nplt.plot(x, likelihood, color = 'blue',label=\"likelihood\")\nplt.plot(x, posterior3, color = 'red',label=\"posterior\")\nplt.legend()"},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"F35C00E89E234C98830A9144987DBD5E","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"可见，合并数据对后验的更新与连续两次对后验更新的结果是一致的。\n\n现在，我们进行一些**练习**：尝试 (1) 不同的数据会对后验产生如何的影响，(2) 通过多次数据对后验分布进行更新"},{"cell_type":"code","execution_count":8,"metadata":{"collapsed":false,"id":"CE134495F00044568F0C94E3F7DEF74E","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"output_type":"error","ename":"TypeError","evalue":"unsupported operand type(s) for -: 'ellipsis' and 'int'","traceback":["\u001b[0;31m---------------------------------------------------------------------------\u001b[0m","\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)","\u001b[0;32m/tmp/ipykernel_172/1333075204.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     15\u001b[0m \u001b[0;31m######################################################\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     16\u001b[0m \u001b[0;31m# 利用已有的函数计算似然\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 17\u001b[0;31m \u001b[0mlikelihood\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mst\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbinom\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpmf\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     18\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     19\u001b[0m \u001b[0;31m# 通过先验与似然，更新后验\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n","\u001b[0;32m/opt/conda/lib/python3.7/site-packages/scipy/stats/_distn_infrastructure.py\u001b[0m in \u001b[0;36mpmf\u001b[0;34m(self, k, *args, **kwds)\u001b[0m\n\u001b[1;32m   3140\u001b[0m         \u001b[0margs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtuple\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmap\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0masarray\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3141\u001b[0m         \u001b[0m_a\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0m_b\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_get_support\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 3142\u001b[0;31m         \u001b[0mk\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0masarray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0mloc\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   3143\u001b[0m         \u001b[0mcond0\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_argcheck\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3144\u001b[0m         \u001b[0mcond1\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mk\u001b[0m \u001b[0;34m>=\u001b[0m \u001b[0m_a\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m&\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mk\u001b[0m \u001b[0;34m<=\u001b[0m \u001b[0m_b\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m&\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_nonzero\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n","\u001b[0;31mTypeError\u001b[0m: unsupported operand type(s) for -: 'ellipsis' and 'int'"]}],"source":"# 加载需要的python包\nimport numpy as np               # numpy 是专门用于数组运算的包\nimport scipy.stats as st         # 从scipy.stats里载入分布函数\nimport matplotlib.pyplot as plt  # matplotlib.pyplot 是专门用于画图的包\n\n# 定义一个Uniform的先验\nx = np.linspace(0,1,100)\nprior = st.uniform.pdf(x, 0, 1)/np.sum(st.uniform.pdf(x,0,1))\n\n######################################################\n# 请填写总抛地球仪的次数\nn = ...\n# 请填写海面朝上的次数\nk = ...\n######################################################\n# 利用已有的函数计算似然\nlikelihood = st.binom.pmf(k,n, x)\n\n# 通过先验与似然，更新后验\nunstd_posterior = prior*likelihood  # 计算分子部分\nposterior = unstd_posterior/np.sum(unstd_posterior) # 计算后验\n\n# 画图\nplt.plot(x, prior , color = 'grey',ls = '--',label=\"prior\")\nplt.plot(x, likelihood, color = 'blue',label=\"likelihood\")\nplt.plot(x, posterior, color = 'red',label=\"posterior\")\nplt.legend()"},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"7A411E7C68CE490E8FD95F269B075FA5","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"(2) 以上一次的后验作为先验，通过新的数据更新后验："},{"cell_type":"code","execution_count":null,"metadata":{"collapsed":false,"id":"D29B9C0E40D24ADFA34F286762AD12A2","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[],"source":"# 以上一次的后验作为先验\nprior = posterior\n\n######################################################\n# 请填写抛地球仪的次数\nn = ...\n# 请填写海面朝上的次数\nk = ...\n######################################################\n# 计算似然\nlikelihood = st.binom.pmf(k,n, x)\n\n# 通过先验与似然，更新后验\nunstd_posterior = prior*likelihood\nposterior = unstd_posterior/np.sum(unstd_posterior)\n\nplt.plot(x, prior , color = 'grey',ls = '--',label=\"prior\")\nplt.plot(x, likelihood, color = 'blue',label=\"likelihood\")\nplt.plot(x, posterior, color = 'red',label=\"posterior\")\nplt.legend()"},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"097491789AF24C2A93EEBC08D87C8B0D","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"最后，我们**练习**设置不同先验对后验的影响。\n- 我们此时假设，海水的覆盖率为0.7或任何0-1的值\n- 并且先验服从正态分布\n- 此时，我们假设收集到的数据为10次数据中有7次海水朝上"},{"cell_type":"code","execution_count":9,"metadata":{"collapsed":false,"id":"2F739B01072649C0B5EB2BE2946EBE58","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"output_type":"error","ename":"TypeError","evalue":"unsupported operand type(s) for -: 'float' and 'ellipsis'","traceback":["\u001b[0;31m---------------------------------------------------------------------------\u001b[0m","\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)","\u001b[0;32m/tmp/ipykernel_172/1844625206.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      7\u001b[0m \u001b[0;31m######################################################\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      8\u001b[0m \u001b[0mx\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinspace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m100\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;31m#生成0到1步长为100的序列\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 9\u001b[0;31m \u001b[0mprior\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mst\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnorm\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpdf\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mmu\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0msd\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mst\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnorm\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpdf\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mmu\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0msd\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     10\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     11\u001b[0m \u001b[0;31m# 计算似然\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n","\u001b[0;32m/opt/conda/lib/python3.7/site-packages/scipy/stats/_distn_infrastructure.py\u001b[0m in \u001b[0;36mpdf\u001b[0;34m(self, x, *args, **kwds)\u001b[0m\n\u001b[1;32m   1868\u001b[0m         \u001b[0margs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtuple\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmap\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0masarray\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1869\u001b[0m         \u001b[0mdtyp\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfind_common_type\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdtype\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfloat64\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1870\u001b[0;31m         \u001b[0mx\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0masarray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mloc\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0mscale\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdtyp\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1871\u001b[0m         \u001b[0mcond0\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_argcheck\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m&\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mscale\u001b[0m \u001b[0;34m>\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1872\u001b[0m         \u001b[0mcond1\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_support_mask\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m&\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mscale\u001b[0m \u001b[0;34m>\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n","\u001b[0;31mTypeError\u001b[0m: unsupported operand type(s) for -: 'float' and 'ellipsis'"]}],"source":"# 定义先验\n######################################################\n# 设置先验的均值：比如0.7\nmu = ...\n# 设置先验的标准差：比如0.1\nsd = ...\n######################################################\nx = np.linspace(0, 1, 100)#生成0到1步长为100的序列\nprior = st.norm.pdf(x,mu,sd)/np.sum(st.norm.pdf(x,mu,sd))\n\n# 计算似然\nlikelihood = st.binom.pmf(7,10, x)\n\n# 通过先验与似然，更新后验\nunstd_posterior = prior*likelihood\nposterior = unstd_posterior/np.sum(unstd_posterior)\n###########绘制先验、似然和后验分布\nplt.plot(x, prior , color = 'grey',ls = '--',label=\"prior\")\nplt.plot(x, likelihood, color = 'blue',label=\"likelihood\")\nplt.plot(x, posterior, color = 'red',label=\"posterior\")\nplt.legend()"},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"6C695C4A9E1A4CFCA5420758CBAEC21B","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"相信大家已经发现了一些**规律**：\n- 后验分布的峰的位置介于似然峰和先验峰之间\n- 随着数据样本量的增加 (抛地球仪次数的增多)，似然函数变得越来越高，越来越窄，对后验的影响越来越大。即后验更接近似然而不是先验。\n\n\n- 不同似然对后验分布的影响\n\t![Image Name](https://cdn.kesci.com/upload/image/rjguv67ra5.png?imageView2/0/w/400/h/400) \n\n- 不同先验对后验分布的影响\n\t![Image Name](https://cdn.kesci.com/upload/image/rjguvhtcgx.png?imageView2/0/w/400/h/400) \n\t\n- 不同数据样本量对后验分布的影响\n\t![Image Name](https://cdn.kesci.com/upload/image/rjguvqbsxi.png?imageView2/0/w/400/h/400)\n</figure>"},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"9DACE3F5DA544733B5EF1F39EC13043D","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"从最后一个图中可以看到，当样本量很大时，后验分布也会很高很窄\n- 这表明我们对于参数 (海水覆盖率) 估计的信心也会更大，即不确定感更少 (uncertarinty，第五节课中关于概率和不确定性的理解)\n\n\n这引出一个问题：当我们获得参数的后验分布时，我们怎么去解释这个分布？\n\n可以认为海水覆盖率为分布的均值？最高点对应的参数值 (信心最高位置对应的参数值)？\n\n我们又怎么结合分布的宽度理解不确定 (uncertarinty)与参数的关系？"},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"03080949F6444530BA6C1071A4E4131F","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"### Point Parameter Estimates"},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"79DD198BDED64CE18A44D843B08BA93A","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"在前面我们知道了如何得到参数的后验分布。\n\n频率论通过样本均值（点估计）进行统计推断 (statistical inference)，但我们该如何通过后验分布进行统计推断?"},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"2EB4A15868A5434E9040A4A5DD0A2696","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"贝叶斯统计中提供了三种对后验分布的点估计量:\n\n- 后验均值 (posterior mean)\n- 后验中位数 (posterior median)\n- 最大后验估计(maximum a posteriori, MAP), 即后验众数 (posterior mode)"},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"5E9EADF64A614B5D8920DD6840EC9E87","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"如果参数后验分布为正态分布，那么后验均值，中位数与众数是接近的。\n- 我们假设参数 (海水覆盖率) 后验分布服从均值为0.7，标值差为0.1的正态分布"},{"cell_type":"code","execution_count":10,"metadata":{"collapsed":false,"id":"9646B9831D064A3DAD38BF4AD996DDBD","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"output_type":"execute_result","data":{"text/plain":"<AxesSubplot:ylabel='Density'>"},"metadata":{},"transient":{}},{"output_type":"display_data","data":{"text/plain":"<Figure size 432x288 with 1 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/9646B9831D064A3DAD38BF4AD996DDBD/rjj8l9tuet.png\">"},"metadata":{"needs_background":"light"},"transient":{}}],"source":"theta_posterior = st.norm.rvs(loc = 0.7, scale = 0.1,size = 50000)  # 生成50000个参数值\n\nsns.kdeplot(theta_posterior) # 使用seaborn的kdeplot绘制参数后验分布"},{"cell_type":"code","execution_count":11,"metadata":{"collapsed":false,"id":"73DCF38D9FFC4CD2BA63179B7C11520E","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"output_type":"stream","name":"stdout","text":"后验均值为： 0.7\n后验中位数为： 0.7\n后验众数为： 0.71\n"}],"source":"import statistics # statistics可以计算数据的平均值，中位数和众数(保留两位小数)\n\nprint(\"后验均值为：\", round(statistics.mean(theta_posterior),2) )\nprint(\"后验中位数为：\", round(statistics.median(theta_posterior),2) )\nprint(\"后验众数为：\",   statistics.mode(theta_posterior.round(2)) )"},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"1568B2667AF0439F9DC245A2B9661568","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"但如果参数的后验分布不是正态分布，而是一个存在偏态和多峰的分布呐？"},{"cell_type":"code","execution_count":12,"metadata":{"collapsed":false,"id":"38F6251FB5754B4190C75D4BA0310003","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"output_type":"execute_result","data":{"text/plain":"<AxesSubplot:ylabel='Density'>"},"metadata":{},"transient":{}},{"output_type":"display_data","data":{"text/plain":"<Figure size 432x288 with 1 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/38F6251FB5754B4190C75D4BA0310003/rjj8lh9l3b.png\">"},"metadata":{"needs_background":"light"},"transient":{}}],"source":"# 假设参数的后验分布服从一个均值为1和5，标准差同为1的多峰分布\ntheta_posterior = np.append(st.norm.rvs(loc = 0.5,scale = 0.05,size = 6000),\nst.norm.rvs(loc = 0.7,scale = 0.05,size = 4000))\n\n# 绘制参数后验分布\nsns.kdeplot(theta_posterior)"},{"cell_type":"code","execution_count":13,"metadata":{"collapsed":false,"id":"FD2ABD7298784E8699D22BBDE76DFFAF","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"output_type":"stream","name":"stdout","text":"后验均值为： 0.58\n后验中位数为： 0.55\n后验众数为： 0.5\n"}],"source":"print(\"后验均值为：\",round(statistics.mean(theta_posterior),2) )\nprint(\"后验中位数为：\",round(statistics.median(theta_posterior),2) )\nprint(\"后验众数为：\", statistics.mode(theta_posterior.round(2)) )"},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"768D04C9A0724F949947A719646BA354","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"其平均值，中位数和众数(MAP)并不相等。\n\n- 其中众数为分布最高点对应的参数值\n- 中位数左右两侧分布的面积各占50%\n- 平均值的位置接近中位数，并且它容易分布形态的影响"},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"C30F28FF2B8A48689A5C16179BD980AE","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"![Image Name](https://cdn.kesci.com/upload/image/rjddhc2wo.png?imageView2/0/w/500/h/500)\n"},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"6533F0539DDB4B658E564473BD0B0202","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"选择哪种点估计进行统计推断？\n\n- 首先后验平均值：均值对分布敏感，更能代表整个分布；并且更贴合测量理论。\n- 次选后验中位数：如果平均值被极值影响，此时最好使用中位数。\n- 慎选后验众数：众数往往存在偏差，并不能很好的代表后验分布。"},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"A49042E926A941569DD0115BDB4B603A","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"### Intervals of Uncertainty"},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"8B6167B9A4DA499C9C0418AC0DEF99A4","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"通过点估计进行统计推断忽略了后验分布的变异性与不确定性 (uncertainty)。\n\n我们通常用区间 (interval)去描述这种不确定性。\n\n- 在贝叶斯的视角中，后验分布的不确定为可信区间(credible interval)。\n- 在频率学派的视角中，通过置信区间 (confidence) 描述数据的不确定性。\n\n然而，我们常常将两者的定义混淆在一起。"},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"1DECC1035BD34C16BE118BD30F6F0DBE","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"#### 贝叶斯可信区间\n\n贝叶斯可信区间 (Bayesian credible intervals)\n- 指参数出现在可信区间内的概率。\n  \n  比如，海水覆盖率95%的可信区间为 [0.6, 0.8m]，那么海水覆盖率为0.7出现在该区间的概率为95%。\n- 描述参数值位置的不确定性。即相对于频率学派认为海水覆盖率是固定值 (比如0.7)，贝叶斯统计可能认为海水覆盖率为0.7的可信度最大。\n- 通过后验分布计算得到的。即相对于频率学派通过抽样数据估计置信区间，可信区间可以通过后验分布进行计算。"},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"2659166B262F453CBD3FE1DC1FA76631","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"贝叶斯置信区间也存在不同的估计指标：\n\n- 后验中心区间(the central posterior interval)\n  - 通过百分位数计算。\n    \n    比如，对于90%后验中心区间，计算后验分布在5%和95%的百分位数即可。\n  - 该区间以参数的中值为中心。\n- 最大后验概率密度区间(the highest density interval, HDI; 或者 the highest posterior density, HPD)\n  - 基于后验概率密度最高占比。\n    \n    比如，90%HDI为后验分布中y轴最高的90%的部分对应的x轴的范围。\n  - 其性质是这个范围中的概率密度永远不低于外部。"},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"4BBC5567E4B14377A9CC22BB3E06E3EA","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"HDI更具代表性，因此通常以HDI作为贝叶斯的可信区间。以下代码展示了如何绘制和计算贝叶斯可信区间。"},{"cell_type":"code","execution_count":14,"metadata":{"collapsed":false,"id":"E4910854B7F547038BE847646938ABCD","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"output_type":"execute_result","data":{"text/plain":"<AxesSubplot:title={'center':'x'}>"},"metadata":{},"transient":{}},{"output_type":"display_data","data":{"text/plain":"<Figure size 432x288 with 1 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/E4910854B7F547038BE847646938ABCD/rjj8lrv9z7.png\">"},"metadata":{"needs_background":"light"},"transient":{}}],"source":"# 假设参数的后验分布服从一个均值为5，标准差为1的正态分布\ntheta_posterior = st.norm.rvs(loc = 5,size = 50000)\n\n# 绘制参数后验分布, 包括其均值与HDI\naz.plot_posterior(theta_posterior,hdi_prob=0.95)"},{"cell_type":"code","execution_count":15,"metadata":{"collapsed":false,"id":"92440FCA0B93434BA3F7C3A0C085E7EE","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"output_type":"stream","name":"stdout","text":"后验分布最大概率密度区间HDI： [3.03172729 6.95505799]\n后验分布的中心区间： [3.0449804  6.97135962]\n"}],"source":"# 通过arviz包的hdi函数计算后验分布的HDI\nprint(\"后验分布最大概率密度区间HDI：\", az.hdi(theta_posterior,hdi_prob=0.95))\n# 通过numpy包的quantile函数计算后验分布的百分位\nprint(\"后验分布的中心区间：\", np.quantile(theta_posterior,[0.025,0.975]))"},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"555CB5A832604FD6BEF5F3DCB73D07FF","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"可见，最大概率密度区间HDI与中心区间是非常接近的。\n\n但是当后验分布为非正态分布时，两种存在区别，比如：\n\n\n![Image Name](https://cdn.kesci.com/upload/image/rjdef48vqj.png?imageView2/0/w/640/h/640)\n\n- 此时，后验中心区间从中位数的两侧开始展开。\n- 而最大概率密度区间HDI从y轴最高的两处开始展开。"},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"D980B999789940528D78019086EE6597","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"#### 我们通过一个练习来了解如何从后验分布中得到点估计和可信区间：\n\n- 首先，我们设置海水的覆盖率的**先验**，比如0-1存在相同的可能。\n- 然后，我们假设收集到的**数据**为10次数据中有7次海水朝上。\n- 接着，我们通过似然和数据更新先验，得到**后验**。\n- 最后，**你的任务是计算后验分布的均值和HDI**。"},{"cell_type":"code","execution_count":16,"metadata":{"collapsed":false,"id":"1259BB5E516941AB93F8EB6A282EC672","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[],"source":"# 定义先验\nx = np.linspace(0,1,100)\nprior = st.uniform.pdf(x,0,1)/np.sum(st.uniform.pdf(x,0,1))\n\n# 计算似然\nlikelihood = st.binom.pmf(7,10, x)\n\n# 通过先验与似然，更新后验\nunstd_posterior = prior*likelihood\nposterior = unstd_posterior/np.sum(unstd_posterior)"},{"cell_type":"code","execution_count":17,"metadata":{"collapsed":false,"id":"62DAFDF2CBFF44378486D3BEA5DB07BA","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"output_type":"error","ename":"TypeError","evalue":"'ellipsis' object is not iterable","traceback":["\u001b[0;31m---------------------------------------------------------------------------\u001b[0m","\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)","\u001b[0;32m/tmp/ipykernel_172/838823818.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      8\u001b[0m \u001b[0;31m# 请使用 statistics.mean() 计算后验分布均值\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      9\u001b[0m \u001b[0;31m###############################################\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 10\u001b[0;31m \u001b[0mposterior_mean\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mstatistics\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmean\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m...\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     11\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"后验均值为：\"\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mround\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mposterior_mean\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n","\u001b[0;32m/opt/conda/lib/python3.7/statistics.py\u001b[0m in \u001b[0;36mmean\u001b[0;34m(data)\u001b[0m\n\u001b[1;32m    304\u001b[0m     \u001b[0mIf\u001b[0m\u001b[0;31m \u001b[0m\u001b[0;31m`\u001b[0m\u001b[0;31m`\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;31m`\u001b[0m\u001b[0;31m`\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0mempty\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mStatisticsError\u001b[0m \u001b[0mwill\u001b[0m \u001b[0mbe\u001b[0m \u001b[0mraised\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    305\u001b[0m     \"\"\"\n\u001b[0;32m--> 306\u001b[0;31m     \u001b[0;32mif\u001b[0m \u001b[0miter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0mdata\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    307\u001b[0m         \u001b[0mdata\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    308\u001b[0m     \u001b[0mn\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n","\u001b[0;31mTypeError\u001b[0m: 'ellipsis' object is not iterable"]}],"source":"# statistics可以计算数据的平均值，中位数和众数(保留两位小数)\nimport statistics\n# arviz 是专门用于贝叶斯统计会绘制的包\nimport arviz as az\n\n###############################################\n# 参数后验分布为 posterior\n# 请使用 statistics.mean() 计算后验分布均值\n###############################################\nposterior_mean = statistics.mean(...)\nprint(\"后验均值为：\",round(posterior_mean,3) )"},{"cell_type":"code","execution_count":18,"metadata":{"collapsed":false,"id":"EC3C9AE4930140C7A03D69335D5A8E75","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"output_type":"error","ename":"ValueError","evalue":"Can only convert xarray dataarray, xarray dataset, dict, netcdf filename, numpy array, pystan fit, pymc3 trace, emcee fit, pyro mcmc fit, numpyro mcmc fit, cmdstan fit csv filename, cmdstanpy fit to InferenceData, not ellipsis","traceback":["\u001b[0;31m---------------------------------------------------------------------------\u001b[0m","\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)","\u001b[0;32m/tmp/ipykernel_172/2541317831.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0;31m###############################################\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0;31m# 通过arviz包的hdi函数计算后验分布的HDI\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 5\u001b[0;31m \u001b[0mposterior_HDI\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0maz\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mhdi\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m...\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mhdi_prob\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0.95\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      6\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"后验分布最大概率密度区间HDI：\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mposterior_HDI\u001b[0m \u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n","\u001b[0;32m/opt/conda/lib/python3.7/site-packages/arviz/stats/stats.py\u001b[0m in \u001b[0;36mhdi\u001b[0;34m(ary, hdi_prob, circular, multimodal, skipna, group, var_names, filter_vars, coords, max_modes, dask_kwargs, **kwargs)\u001b[0m\n\u001b[1;32m    461\u001b[0m         \u001b[0mary\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexpand_dims\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mary\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    462\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 463\u001b[0;31m     \u001b[0mary\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mconvert_to_dataset\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mary\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mgroup\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mgroup\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    464\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mcoords\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    465\u001b[0m         \u001b[0mary\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mget_coords\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mary\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcoords\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n","\u001b[0;32m/opt/conda/lib/python3.7/site-packages/arviz/data/converters.py\u001b[0m in \u001b[0;36mconvert_to_dataset\u001b[0;34m(obj, group, coords, dims)\u001b[0m\n\u001b[1;32m    177\u001b[0m     \u001b[0mxarray\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mDataset\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    178\u001b[0m     \"\"\"\n\u001b[0;32m--> 179\u001b[0;31m     \u001b[0minference_data\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mconvert_to_inference_data\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobj\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mgroup\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mgroup\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcoords\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcoords\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdims\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdims\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    180\u001b[0m     \u001b[0mdataset\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mgetattr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0minference_data\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mgroup\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    181\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mdataset\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n","\u001b[0;32m/opt/conda/lib/python3.7/site-packages/arviz/data/converters.py\u001b[0m in \u001b[0;36mconvert_to_inference_data\u001b[0;34m(obj, group, coords, dims, **kwargs)\u001b[0m\n\u001b[1;32m    131\u001b[0m         raise ValueError(\n\u001b[1;32m    132\u001b[0m             \"Can only convert {} to InferenceData, not {}\".format(\n\u001b[0;32m--> 133\u001b[0;31m                 \u001b[0;34m\", \"\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mjoin\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mallowable_types\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mobj\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__class__\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__name__\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    134\u001b[0m             )\n\u001b[1;32m    135\u001b[0m         )\n","\u001b[0;31mValueError\u001b[0m: Can only convert xarray dataarray, xarray dataset, dict, netcdf filename, numpy array, pystan fit, pymc3 trace, emcee fit, pyro mcmc fit, numpyro mcmc fit, cmdstan fit csv filename, cmdstanpy fit to InferenceData, not ellipsis"]}],"source":"###############################################\n# 请使用 az.hdi() 计算后验分布HDI\n###############################################\n# 通过arviz包的hdi函数计算后验分布的HDI\nposterior_HDI = az.hdi(..., hdi_prob=0.95)\nprint(\"后验分布最大概率密度区间HDI：\", posterior_HDI )"},{"cell_type":"code","execution_count":19,"metadata":{"collapsed":false,"id":"B4EF3329E7534F0A8783291B319BB0F9","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"output_type":"error","ename":"ValueError","evalue":"Can only convert xarray dataarray, xarray dataset, dict, netcdf filename, numpy array, pystan fit, pymc3 trace, emcee fit, pyro mcmc fit, numpyro mcmc fit, cmdstan fit csv filename, cmdstanpy fit to InferenceData, not ellipsis","traceback":["\u001b[0;31m---------------------------------------------------------------------------\u001b[0m","\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)","\u001b[0;32m/tmp/ipykernel_172/3182002862.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      2\u001b[0m \u001b[0;31m# 请使用 az.plot_posterior()函数绘制后验分布图像\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0;31m###############################################\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 4\u001b[0;31m \u001b[0maz\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot_posterior\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m...\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m","\u001b[0;32m/opt/conda/lib/python3.7/site-packages/arviz/plots/posteriorplot.py\u001b[0m in \u001b[0;36mplot_posterior\u001b[0;34m(data, var_names, filter_vars, transform, coords, grid, figsize, textsize, hdi_prob, multimodal, skipna, round_to, point_estimate, group, rope, ref_val, kind, bw, circular, bins, ax, backend, backend_kwargs, show, **kwargs)\u001b[0m\n\u001b[1;32m    202\u001b[0m         \u001b[0;34m>>\u001b[0m\u001b[0;34m>\u001b[0m \u001b[0maz\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot_posterior\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvar_names\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'mu'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mhdi_prob\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m.75\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    203\u001b[0m     \"\"\"\n\u001b[0;32m--> 204\u001b[0;31m     \u001b[0mdata\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mconvert_to_dataset\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mgroup\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mgroup\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    205\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mtransform\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    206\u001b[0m         \u001b[0mdata\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtransform\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n","\u001b[0;32m/opt/conda/lib/python3.7/site-packages/arviz/data/converters.py\u001b[0m in \u001b[0;36mconvert_to_dataset\u001b[0;34m(obj, group, coords, dims)\u001b[0m\n\u001b[1;32m    177\u001b[0m     \u001b[0mxarray\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mDataset\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    178\u001b[0m     \"\"\"\n\u001b[0;32m--> 179\u001b[0;31m     \u001b[0minference_data\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mconvert_to_inference_data\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobj\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mgroup\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mgroup\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcoords\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcoords\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdims\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdims\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    180\u001b[0m     \u001b[0mdataset\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mgetattr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0minference_data\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mgroup\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    181\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mdataset\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n","\u001b[0;32m/opt/conda/lib/python3.7/site-packages/arviz/data/converters.py\u001b[0m in \u001b[0;36mconvert_to_inference_data\u001b[0;34m(obj, group, coords, dims, **kwargs)\u001b[0m\n\u001b[1;32m    131\u001b[0m         raise ValueError(\n\u001b[1;32m    132\u001b[0m             \"Can only convert {} to InferenceData, not {}\".format(\n\u001b[0;32m--> 133\u001b[0;31m                 \u001b[0;34m\", \"\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mjoin\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mallowable_types\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mobj\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__class__\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__name__\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    134\u001b[0m             )\n\u001b[1;32m    135\u001b[0m         )\n","\u001b[0;31mValueError\u001b[0m: Can only convert xarray dataarray, xarray dataset, dict, netcdf filename, numpy array, pystan fit, pymc3 trace, emcee fit, pyro mcmc fit, numpyro mcmc fit, cmdstan fit csv filename, cmdstanpy fit to InferenceData, not ellipsis"]}],"source":"###############################################\n# 请使用 az.plot_posterior()函数绘制后验分布图像\n###############################################\naz.plot_posterior(...)"},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"22A3CB1362554E5C95051E1EBE1E5827","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"#### 频率学置信区间 vs. 贝叶斯可信区间\n\n在频率学派统计中，比如我们抛地球仪100次，其中有70次海面朝上，其95%的置信区间为[0.6, 0.8]\n\n此时我们能认为，海水覆盖率=0.7在这个区间的概率为95%吗, 或者海水覆盖率在[0.6, 0.8]之外时的概率为5%？"},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"AE74BD41E9B14BF499A0721B13C8B045","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"上面的理解仿佛是显而易见的，但事实却是：\n\n- 错误的理解：真实海水覆盖率在这个区间的概率为95%。\n- 正确的解释：当我们收集无数次数据(即获得无数个抛100次地球仪的数据)，然后根据数据计算的无数个置信区间，这无数个置信区间包括真实海水覆盖率=0.7在这个区间的概率为95%的概率为95%。\n\n显然，频率学派对于置信区间的定义是难以理解的。\n\n而那个错误的理解更容易让人接受，这个错误的解释恰好解释贝叶斯可信区间的含义。"},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"E8D5D232BDB547F0AD1A3AFF079997C4","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"在现实中，我们只从总体中抽取一个样本，这使得我们没有办法知道置信区间是否包含真正的参数值。\n\n这意味着我们只能通过置信区间否定正值可能不在这个区间内（零假设检验）, 但是我们无法得到真值在这个范围的概率。\n\n并且基于单个样本推测的置信区间易受到极端值或者抽样误差的影响。\n\n置信区间表示我们从数据推测真值的区间的不确定性，而不是有关感兴趣参数出现的概率。"},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"979B5DBC2E6C43B6BB36B38A58DB6F25","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"置信区间与可信区间反应了对于不确定的不同看法：\n- 置信区间认为参数是固定不变的真值，而不确定的是置信区间本身，因此置信区间描述的是该样本推测出特定置信区间的可能性。\n- 可信区间认为观测数据是不变的，而参数是变化的，因此不确定性描述的是参数变化的可能性。"},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"A0D928073D44404390D6AE28B3BABC6D","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"与频率学派的置信区间相比：\n\n可信区间更容易理解，即可信区间就是指参数或者现象出现的范围。\n\n其不确定性通过概率表示，代表了我们对于参数出现的信心。\n\n而不是如频率学派的“在无数次抽样后，得到的区间包括真值的范围和概率”。"},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"766BB64DE8CE4CC1999FF1DA7B2D8203","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"## Part 3: From Posterior to Predictions by Sampling"},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"368270DD780840BFA12605FC812D4E98","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"后验分布的另一大优势是可以进行预测 (prediction). \n\n具体步骤\n- 我们可以从参数后验分布中选取部分参数。\n  \n  $θ_{i} \\sim p(θ | data)$, 其中$p(θ | data)$为参数后验分布，$θ_{i}$为从后验分布中选取的i个参数。\n- 通过参数生成数据。\n  \n  $data_{i} \\sim p(data | θ_{i})$, 其中$p(data | θ_{i})$为似然函数 (likelihood)，$data_{i}$为从生成的第i个数据。\n\n- 通过重复以上两个步骤，我们可以得到生成数据的分布，即后验预测分布 (posterior predictive distribution)。"},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"A0435B4EFBB149F2A1CD92198FD72600","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"后验预测分布（预测数据）与真实数据的差异体现在**两种变异** (uncertainty)上：\n1. 参数抽取的变异。比如，我们抽取的参数可能来自参数后验分布的尾部。\n2. 数据生成的变异。比如，我们生成的上述可能来自似然函数的尾部。"},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"FA0E9E393B1D476C8A855AFFFD224AF2","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"#### 接下来我们通过代码展示如何通过后验参数生成数据。\n\n以探究地球海水覆盖率为例。\n- 假如我们得到了关于地球海水覆盖率的后验概率分布。\n- 其中，海水覆盖率为0.7，并且服从均值为0.7，标准差为0.1的正态分布。"},{"cell_type":"code","execution_count":20,"metadata":{"collapsed":false,"id":"2101DE8BD165421CB3AFB7CD84E132DE","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"output_type":"execute_result","data":{"text/plain":"Text(0.5, 0, 'Sea water coverage')"},"metadata":{},"transient":{}},{"output_type":"display_data","data":{"text/plain":"<Figure size 432x288 with 1 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/2101DE8BD165421CB3AFB7CD84E132DE/rjj8mp5ogb.png\">"},"metadata":{"needs_background":"light"},"transient":{}}],"source":"# 从scipy.stats里载入二项分布函数\nfrom scipy.stats import norm\n# matplotlib.pyplot 是专门用于画图的包\nimport matplotlib.pyplot as plt\n# seaborn是专门用于绘制统计数据的包\nimport seaborn as sns\n\n# 生成参数后验分布\ntheta_posterior = norm.rvs(loc = 0.7, scale=0.1,size = 50000)\n\n# 绘制分布图\nsns.kdeplot(theta_posterior)\nplt.vlines(0.7,ymin=0,ymax=4,color=\"red\")\nplt.xlabel(\"Sea water coverage\")"},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"C2CF445F5337478CBE0B0371005B6AC8","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"执行步骤一：从后验分布中抽取采样"},{"cell_type":"code","execution_count":21,"metadata":{"collapsed":false,"id":"B87BC97CDCA44E9F8DD78A191096A334","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"output_type":"stream","name":"stdout","text":"抽取的参数值为： [0.78521256]\n"}],"source":"# numpy 是专门用于数组运算的包\nimport numpy as np\n\n# random.choice函数可以实现从后验分布中抽取参数，其中size为抽取的参数的数量\ntheta_i = np.random.choice(theta_posterior, size = 1)\nprint(\"抽取的参数值为：\",theta_i)"},{"cell_type":"code","execution_count":22,"metadata":{"collapsed":false,"id":"672BE514929540EEAEEF955222AC0E76","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"output_type":"error","ename":"TypeError","evalue":"int() argument must be a string, a bytes-like object or a number, not 'ellipsis'","traceback":["\u001b[0;31m---------------------------------------------------------------------------\u001b[0m","\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)","\u001b[0;32m/tmp/ipykernel_172/1734117540.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     15\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     16\u001b[0m \u001b[0;31m# random.choice函数可以实现从后验分布中抽取参数，其中size为抽取的参数的数量\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 17\u001b[0;31m \u001b[0mtheta_is\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrandom\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mchoice\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtheta_posterior\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0msize\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mnumber\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     18\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"抽取的参数值为：\"\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mtheta_is\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n","\u001b[0;32mmtrand.pyx\u001b[0m in \u001b[0;36mnumpy.random.mtrand.RandomState.choice\u001b[0;34m()\u001b[0m\n","\u001b[0;32m<__array_function__ internals>\u001b[0m in \u001b[0;36mprod\u001b[0;34m(*args, **kwargs)\u001b[0m\n","\u001b[0;32m/opt/conda/lib/python3.7/site-packages/numpy/core/fromnumeric.py\u001b[0m in \u001b[0;36mprod\u001b[0;34m(a, axis, dtype, out, keepdims, initial, where)\u001b[0m\n\u001b[1;32m   3029\u001b[0m     \"\"\"\n\u001b[1;32m   3030\u001b[0m     return _wrapreduction(a, np.multiply, 'prod', axis, dtype, out,\n\u001b[0;32m-> 3031\u001b[0;31m                           keepdims=keepdims, initial=initial, where=where)\n\u001b[0m\u001b[1;32m   3032\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3033\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n","\u001b[0;32m/opt/conda/lib/python3.7/site-packages/numpy/core/fromnumeric.py\u001b[0m in \u001b[0;36m_wrapreduction\u001b[0;34m(obj, ufunc, method, axis, dtype, out, **kwargs)\u001b[0m\n\u001b[1;32m     85\u001b[0m                 \u001b[0;32mreturn\u001b[0m \u001b[0mreduction\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0maxis\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0maxis\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mout\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mout\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mpasskwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     86\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 87\u001b[0;31m     \u001b[0;32mreturn\u001b[0m \u001b[0mufunc\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreduce\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobj\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mout\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mpasskwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     88\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     89\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n","\u001b[0;31mTypeError\u001b[0m: int() argument must be a string, a bytes-like object or a number, not 'ellipsis'"]}],"source":"######################################################\n# 练习阶段\n#\n# 目的：了解后验参数抽取过程中的变异\n# \n# Tips: 通过抽取多个参数值观察，是否抽取的参数值存在围绕均值0.7的变异\n######################################################\n# numpy 是专门用于数组运算的包\nimport numpy as np\n\n######################################################\n# 请填写抽取参数的数量\nnumber = ...\n######################################################\n\n# random.choice函数可以实现从后验分布中抽取参数，其中size为抽取的参数的数量\ntheta_is = np.random.choice(theta_posterior,size=number)\nprint(\"抽取的参数值为：\",theta_is)"},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"BCD5646FB1FD48E0B3002FCCECEFBA00","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"执行步骤二：通过参数值生成数据。\n\n假设我们抛了100次地球仪，加上我们已经知道地球海面的覆盖率（即步骤一抽到的参数）。\n\n我们可以利用上面的参数去预测，100次数据中有多少次海面朝上的数据。"},{"cell_type":"code","execution_count":23,"metadata":{"collapsed":false,"id":"F6AB4879CB84409FA932FE438FB17D96","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"output_type":"stream","name":"stdout","text":"海面朝上的次数： 74\n"}],"source":"# 从scipy.stats里载入二项分布函数\nfrom scipy.stats import binom,norm\n\n# 假设抛100次地球仪\nn = 100\n# 带入抽取的参数，生成假数据\ndata_i = binom.rvs(n,theta_i)\nprint(\"海面朝上的次数：\", data_i)"},{"cell_type":"code","execution_count":24,"metadata":{"collapsed":false,"id":"1CA1539A0DEF4B599848AC4FE4218B5E","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":true,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"output_type":"error","ename":"TypeError","evalue":"'ellipsis' object cannot be interpreted as an integer","traceback":["\u001b[0;31m---------------------------------------------------------------------------\u001b[0m","\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)","\u001b[0;32m/tmp/ipykernel_172/2682088426.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     17\u001b[0m \u001b[0mn\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m100\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     18\u001b[0m \u001b[0;31m# 带入抽取的参数，生成假数据\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 19\u001b[0;31m \u001b[0mdata_is\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mbinom\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrvs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mtheta_i\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0msize\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mnumber\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     20\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"海面朝上的次数：\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdata_is\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n","\u001b[0;32m/opt/conda/lib/python3.7/site-packages/scipy/stats/_distn_infrastructure.py\u001b[0m in \u001b[0;36mrvs\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m   3115\u001b[0m         \"\"\"\n\u001b[1;32m   3116\u001b[0m         \u001b[0mkwargs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'discrete'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 3117\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0msuper\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrvs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   3118\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3119\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mpmf\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mk\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwds\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n","\u001b[0;32m/opt/conda/lib/python3.7/site-packages/scipy/stats/_distn_infrastructure.py\u001b[0m in \u001b[0;36mrvs\u001b[0;34m(self, *args, **kwds)\u001b[0m\n\u001b[1;32m   1090\u001b[0m             \u001b[0mvals\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_rvs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1091\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1092\u001b[0;31m             \u001b[0mvals\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_rvs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msize\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0msize\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mrandom_state\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mrandom_state\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1093\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1094\u001b[0m         \u001b[0mvals\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mvals\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mscale\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mloc\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n","\u001b[0;32m/opt/conda/lib/python3.7/site-packages/scipy/stats/_discrete_distns.py\u001b[0m in \u001b[0;36m_rvs\u001b[0;34m(self, n, p, size, random_state)\u001b[0m\n\u001b[1;32m     50\u001b[0m     \"\"\"\n\u001b[1;32m     51\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m_rvs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mn\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mp\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msize\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mrandom_state\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 52\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0mrandom_state\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbinomial\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mp\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msize\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     53\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     54\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m_argcheck\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mn\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mp\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n","\u001b[0;32mmtrand.pyx\u001b[0m in \u001b[0;36mnumpy.random.mtrand.RandomState.binomial\u001b[0;34m()\u001b[0m\n","\u001b[0;31mTypeError\u001b[0m: 'ellipsis' object cannot be interpreted as an integer"]}],"source":"######################################################\n# 练习阶段\n#\n# 目的：了解数据生成过程中的变异\n# \n# Tips: 通过生成多个数据观察，是否生成的数据存在某些规律\n######################################################\n# 从scipy.stats里载入二项分布函数\nfrom scipy.stats import binom\n\n######################################################\n# 请填写抽取参数的数量\nnumber = ...\n######################################################\n\n# 假设抛100次地球仪\nn = 100\n# 带入抽取的参数，生成假数据\ndata_is = binom.rvs(n,theta_i,size=number)\nprint(\"海面朝上的次数：\", data_is)"},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"233AB3B76369478C9D83B48B4FBBD0BA","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"至此，我们已经了解了如何通过参数后验分布预测对数据进行预测。并且了解了数据生成过程中的两种变异。"},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"6ACA68332757409BB51668DAC7EF6A03","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"#### 最后，我们通过绘图的方式展示后验预测分布与真实数据的关系。\n\n首先，我们假设我们获得的真实数据 (observed_data) 为：抛100地球仪获得了70次海面朝上。\n\n然后，假设我们已经根据先验和数据得到地球海水覆盖率的后验概率 (theta_posterior)。\n\n接着，我们通过后验概率分布，生成预测数据 (predictive_data)。\n\n最后，我们比较**预测数据**与**真实数据**的差异。"},{"cell_type":"code","execution_count":25,"metadata":{"collapsed":false,"id":"913DFFC4A3C24E92B7B68CDFDE2894D4","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[],"source":"# numpy 是专门用于数组运算的包\nimport numpy as np\n# 从scipy.stats里载入二项分布函数\nfrom scipy.stats import binom,norm\n\n# random.choice函数可以实现从后验分布中抽取参\n# 我们从后验分布中抽取50个参数\ntheta_is = np.random.choice(theta_posterior,size=50)\n\n# 假设抛100次地球仪\nn = 100\n# 带入抽取的参数，生成假数据\npredictive_data = np.array([])\nfor theta_i in theta_is:\n  # 每次生成100份数据\n  predictive_data = np.append(predictive_data,binom.rvs(n,theta_i,size=100))"},{"cell_type":"code","execution_count":26,"metadata":{"collapsed":false,"id":"26C4D58C03204692909FD49C7F123AD2","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"output_type":"execute_result","data":{"text/plain":"Text(0.5, 0, 'The number of times the sea is facing up')"},"metadata":{},"transient":{}},{"output_type":"display_data","data":{"text/plain":"<Figure size 432x288 with 1 Axes>","text/html":"<img src=\"https://cdn.kesci.com/upload/rt/26C4D58C03204692909FD49C7F123AD2/rjj8ndiucl.png\">"},"metadata":{"needs_background":"light"},"transient":{}}],"source":"sns.kdeplot(predictive_data)  # 绘制预测数据的概率密度图\n\n# 对比真实数据与预测数据\nobserved_data = 70\npredictive_data_mean = np.mean(predictive_data) \nplt.vlines(observed_data,ymin=0,ymax=0.038,color=\"red\")          # 红色竖线为数据真实值\nplt.vlines(predictive_data_mean,ymin=0,ymax=0.038,color=\"black\") # 黑色竖线为预测数据的平均值\nplt.xlabel(\"The number of times the sea is facing up\")"},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"960A5497DB7346FB9049F7C0C7D7AA47","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"可以看到，数据预测的均值（黑线）与真实数据（红线）非常接近。"},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"FFEB92CAF6114003985AD74345045CD9","jupyter":{},"notebookId":"6343d1cebdb14aa731d093af","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"总结：\n- 了解了贝叶斯公式中分母（边际似然）对于后验的影响\n- 通过参数后验分布进行统计推断，包括通过点估计的三种方式与可信区间的两种方式。\n- 比较了频率学派置信区间与贝叶斯可信区间的区别\n- 通过代码练习如何从后验分布中生成预测数据"}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python","nbconvert_exporter":"python","file_extension":".py","version":"3.5.2","pygments_lexer":"ipython3"}},"nbformat":4,"nbformat_minor":2}